#Set up

rm(list=ls())

#Set your working directory - this will be the starting folder where all other folders are created
setwd ("C:/Users/cecidy/Google Drive/ZH_analysis/")

source("ZHfunctions.R")
options(scipen=9999)

######################################################################################################

#Analysis for the paper:

#"Assessing multi-dimensional sustainability: lessons from Brazil's social protection programs"

#Robustness tests


#Robustness test 1: Spatial autocorrelation
#Robustness test 2: Endogeneity

######################################################################################################
######################################################################################################

#Functions sourced from ZHfunctions.R, in combination with mapply are used to carry out the analysis

#There are 24 programme-outcome combinations and thus 24 models
#Arguments needed for the functions are read from the lists below
#The lists are sourced again before each function is run
#Analysis starts around line 250

######################################################################################################
######################################################################################################

#starttest

#A list of directory where the script for each model lie
#If all scripts lie in the same folder, repeat 24 times (one per programme-outcome combination)

Folderlist <- list("ProgrammeOutcome/ZH/Kcal/",
                   "ProgrammeOutcome/BF/Kcal/",
                   "ProgrammeOutcome/PRONAF/Kcal/",
                   
                   "ProgrammeOutcome/ZH/Protein/",
                   "ProgrammeOutcome/BF/Protein/",
                   "ProgrammeOutcome/PRONAF/Protein/",
                   
                   "ProgrammeOutcome/ZH/MPI00to10/",
                   "ProgrammeOutcome/BF/MPI00to10/",
                   "ProgrammeOutcome/PRONAF/MPI00to10/",
                   
                   "ProgrammeOutcome/ZH/MPI04to13/",
                   "ProgrammeOutcome/BF/MPI04to13/",
                   "ProgrammeOutcome/PRONAF/MPI04to13/",
                   
                   "ProgrammeOutcome/ZH/Childmal/",
                   "ProgrammeOutcome/BF/Childmal/",
                   "ProgrammeOutcome/PRONAF/Childmal/",
                   
                   "ProgrammeOutcome/ZH/InfMort00to10/",
                   "ProgrammeOutcome/BF/InfMort00to10/",
                   "ProgrammeOutcome/PRONAF/InfMort00to10/",
                   
                   "ProgrammeOutcome/ZH/InfMort04to13/",
                   "ProgrammeOutcome/BF/InfMort04to13/",
                   "ProgrammeOutcome/PRONAF/InfMort04to13/",
                   
                   "ProgrammeOutcome/ZH/NatVeg/",
                   "ProgrammeOutcome/BF/NatVeg/",
                   "ProgrammeOutcome/PRONAF/NatVeg/")

Scriptlist <- list("1.ZHandKcal_core.R",
                   "2.BFandKcal_core.R",
                   "3.PRONAFandKcal_core.R",
                   "4.ZHandProtein_core.R",
                   "5.BFandProtein_core.R",
                   "6.PRONAFandProtein_core.R",
                   "7.ZHandMPI00to10_core.R",
                   "8.BFandMPI00to10_core.R",
                   "9.PRONAFandMPI00to10_core.R",
                   "10.ZHandMPI04_core.R",
                   "11.BFandMPI04_core.R",
                   "12.PRONAFandMPI04_core.R",
                   "13.ZHandChildmal_core.R",
                   "14.BFandChildmal_core.R",
                   "15.PRONAFandChildmal_core.R",
                   "16.ZHandInfmort00to10_core.R",
                   "17.BFandInfmort00to10_core.R",
                   "18.PRONAFandInfmort00to10_core.R",
                   "19.ZHandInfmort04to13_core.R",
                   "20.BFandInfmort04to13_core.R",
                   "21.PRONAFandInfmort04to13_core.R",
                   "22.ZHandNatVeg_core.R",
                   "23.BFandNatVeg_core.R",
                   "24.PRONAFandNatVeg_core.R")

#Bookmarks in scripts to indicate rows from where to start to run script
startList <- list("#start2_","#start2_","#start2_",
                  "#start2_","#start2_","#start2_",
                  "#start2_","#start2_","#start2_",
                  "#start2_","#start2_","#start2_",
                  "#temp1","#temp1","#temp1",
                  "#temp1","#temp1","#temp1",
                  "#temp1","#temp1","#temp1",
                  "#start2_","#start2_","#start2_")

#Bookmarks in scripts to indicate rows from where to stop to run script
endList <- list("#fin2_","#fin2_","#fin2_",
                "#fin2_","#fin2_","#fin2_",
                "#fin2_","#fin2_","#fin2_",
                "#fin2_","#fin2_","#fin2_",
                "#temp2","#temp2","#temp2",
                "#temp2","#temp2","#temp2",
                "#temp2","#temp2","#temp2",
                "#fin2_","#fin2_","#fin2_")

#YearSampList, TimeperiodList and ProgList are lists used to read in the correct
#election variable used in each model.

YearSampList <- list("2004","2004","2004","2004","2004","2004",
                     "2000","2000","2000","2004","2004","2004",
                     "2004","2004","2004", 
                     "2000","2000","2000","2004","2004","2004",
                     "2004","2004","2004")


TimeperiodList <- list("2004to2013","2004to2013","2004to2013","2004to2013","2004to2013","2004to2013",
                       "2000to2010","2000to2010","2000to2010","2004to2013","2004to2013","2004to2013",
                       "2004to2013","2004to2013","2004to2013",
                       "2000to2010","2000to2010","2000to2010","2004to2013","2004to2013","2004to2013",
                       "2004to2013","2004to2013","2004to2013")

ProgList <- list("ZH","BF","PRONAF","ZH","BF","PRONAF",
                 "ZH","BF","PRONAF","ZH","BF","PRONAF",
                 "ZH","BF","PRONAF",
                 "ZH","BF","PRONAF","ZH","BF","PRONAF",
                 "ZH","BF","PRONAF")

#fitnamesave - this will be the same for ALL!

#List of folders where the outcomes (e.g. CBGPS weights) should be saved
#These are put in separate folder per programme-outcome, and named the same, so that the files can easily be read in again
#and used in another step
FolderSavelist <- list("ProgrammeOutcome/ZH/Kcal/",
                       "ProgrammeOutcome/BF/Kcal/",
                       "ProgrammeOutcome/PRONAF/Kcal/",
                       
                       "ProgrammeOutcome/ZH/Protein/",
                       "ProgrammeOutcome/BF/Protein/",
                       "ProgrammeOutcome/PRONAF/Protein/",
                       
                       "ProgrammeOutcome/ZH/MPI00to10/",
                       "ProgrammeOutcome/BF/MPI00to10/",
                       "ProgrammeOutcome/PRONAF/MPI00to10/",
                       
                       "ProgrammeOutcome/ZH/MPI04to13/",
                       "ProgrammeOutcome/BF/MPI04to13/",
                       "ProgrammeOutcome/PRONAF/MPI04to13/",
                       
                       "ProgrammeOutcome/ZH/Childmal/",
                       "ProgrammeOutcome/BF/Childmal/",
                       "ProgrammeOutcome/PRONAF/Childmal/",
                       
                       "ProgrammeOutcome/ZH/InfMort00to10/",
                       "ProgrammeOutcome/BF/InfMort00to10/",
                       "ProgrammeOutcome/PRONAF/InfMort00to10/",
                       
                       "ProgrammeOutcome/ZH/InfMort04to13/",
                       "ProgrammeOutcome/BF/InfMort04to13/",
                       "ProgrammeOutcome/PRONAF/InfMort04to13/",
                       
                       "ProgrammeOutcome/ZH/NatVeg/",
                       "ProgrammeOutcome/BF/NatVeg/",
                       "ProgrammeOutcome/PRONAF/NatVeg/")

#Three different variations of biome are used, depending on the model,
#Biomefor2004Analysis is for models where municipality borders are adjusted to 2004 boundaries, i.e. 2004-2013 time-frame models
#Biomefor2000Analysis is for models where municipality borders are adjusted to 2000 boundaries, i.e. 2000-2010 time-frame models
#Biomefor2004AnalysisNoHybrid is used in Natural vegetation models, are adjusted to 2004 boundaris, and only have main biome categories (no transition ones)

biomeL <- list("Biomefor2004Analysis","Biomefor2004Analysis","Biomefor2004Analysis",
               "Biomefor2004Analysis","Biomefor2004Analysis","Biomefor2004Analysis",
               "Biomefor2000Analysis","Biomefor2000Analysis","Biomefor2000Analysis",
               "Biomefor2004Analysis","Biomefor2004Analysis","Biomefor2004Analysis",
               
               "Biomefor2004Analysis","Biomefor2004Analysis","Biomefor2004Analysis",
               "Biomefor2000Analysis","Biomefor2000Analysis","Biomefor2000Analysis",
               "Biomefor2004Analysis","Biomefor2004Analysis","Biomefor2004Analysis",
               
               "Biomefor2004AnalysisNoHybrid","Biomefor2004AnalysisNoHybrid","Biomefor2004AnalysisNoHybrid")


#The item number for the best fit models used in analysis (to run robustness tests on)
NrList <- list(3,3,3,
               3,3,3,
               
               3,3,3,
               1,1,3,
               
               4,4,4,
               1,1,1,
               4,4,4,
               
               3,3,3)

#The Brazil municipality border shapefile used to test for spatial autocorrelation
#MU2001_proj are borders used for the time frame 2004-2013, while MU2000_proj are borders used for the time frame 2000-2010 
ShapeFile <- list("MU2001_proj","MU2001_proj","MU2001_proj",
                  "MU2001_proj","MU2001_proj","MU2001_proj",
                  
                  "MU2000_proj","MU2000_proj","MU2000_proj",
                  "MU2001_proj","MU2001_proj","MU2001_proj",
                  
                  "MU2001_proj","MU2001_proj","MU2001_proj",
                  "MU2000_proj","MU2000_proj","MU2000_proj",
                  "MU2001_proj","MU2001_proj","MU2001_proj",
                  
                  "MU2001_proj","MU2001_proj","MU2001_proj")

#Select whether need to exclude influential points (sqrtFunc) or not (sqrtFuncEmpty) for the model
sqrtL <- list(sqrtFuncEmpty,sqrtFuncEmpty,sqrtFuncEmpty,
              sqrtFuncEmpty,sqrtFuncEmpty,sqrtFuncEmpty,
              sqrtFuncEmpty,sqrtFuncEmpty,sqrtFuncEmpty,
              sqrtFuncEmpty,sqrtFuncEmpty,sqrtFuncEmpty,
              
              sqrtFunc,sqrtFunc,sqrtFunc,
              sqrtFunc,sqrtFunc,sqrtFunc,
              sqrtFunc,sqrtFunc,sqrtFunc,
              
              sqrtFuncEmpty,sqrtFuncEmpty,sqrtFuncEmpty)

#Only need number and model for the QP/NB ones - these are the ones that the influential points are removed from
#only pronaf-childmal was state linear in the non-robust model
nrL <- list(NA,NA,NA,
            NA,NA,NA,
            NA,NA,NA,
            NA,NA,NA,
            
            4,4,6,
            1,1,1,
            4,4,4,
            
            NA,NA,NA)

#endtest

###############################################################################################################

#1. Test for spatial autocorrelation

##############################################################################################################

#Test for spatial autocorrelation for one model

test <- NewestMoranALL2019Elec("ProgrammeOutcome/ZH/Kcal/","1.ZHandKcal_core.R",
                               "#start2_","#fin2_","2004","2004to2013","ZH",sqrtFuncEmpty,"ProgrammeOutcome/ZH/Kcal/",
                               "AllCoreMods_W_Election_Full.rda",NA,"Biomefor2004Analysis","AllCoreMods_W_Election_NoInfl.rda",3,"MU2001_proj")

#if get warning can ignore, its due to observations missing 
#from online: zero sum general weights is when the general weights sum is zero for 
#this region. It's a warning, not an error - you've set zero.policy to 
#TRUE, so the output is OK assuming that you accept no-neighbour observations.

#######################################################################################################

#For all core models
sourcePartial("ZH_RobustnessTests.R",startTag="#starttest",endTag="#endtest")

MoranAutoRes1 <- mapply(function(x,y,z,w,a,b,c,d,e,f,g,h,i) NewestMoranALL2019Elec(x,y,z,w,a,b,c,d,e,"AllCoreMods_W_Election_Full.rda",f,g,
                                                                                   "AllCoreMods_W_Election_NoInfl.rda",h,i),Folderlist,Scriptlist,
                        startList,endList,YearSampList,TimeperiodList,ProgList,sqrtL,FolderSavelist,nrL,biomeL,NrList,ShapeFile,
                        SIMPLIFY = FALSE)

##########################################################################################################################

#Prepare files so that I can just read all of them back in (from .csv)
sourcePartial("ZH_RobustnessTests.R",startTag="#starttest",endTag="#endtest")


#CORE

#1. Propensity model
Together <- paste(FolderSavelist,"AllCoreMods_W_Election_NoInfl.rda",sep="",paste("_",NrList,sep="",
                                                                                  paste("_","MoranTablePropensityWElection",sep="",paste(".csv",sep=""))))

PropensityTogetherCSV <- lapply(Together, function(y) read.csv(y, header=TRUE,sep=","))
PropensityAll <- data.frame(do.call("rbind",PropensityTogetherCSV))

TogetherName <- paste(FolderSavelist,"AllCoreMods_W_Election_NoInfl.rda",sep="",paste("_",NrList,sep="",paste(".csv",sep="")))
TogetherName <- rep(TogetherName,each=2)

Propensity <- "PropensityMod"
PropensityAllFin <- cbind(TogetherName,Propensity,PropensityAll)

#2. Outcome model

Together <- paste(FolderSavelist,"AllCoreMods_W_Election_NoInfl.rda",sep="",paste("_",NrList,sep="",paste("MoranTableOutcomeWElection",sep="",paste(".csv",sep=""))))

OutcomeTogetherCSV <- lapply(Together, function(y) read.csv(y, header=TRUE,sep=","))
OutcomeAll <- data.frame(do.call("rbind",OutcomeTogetherCSV))

TogetherName <- paste(FolderSavelist,"AllCoreMods_W_Election_NoInfl.rda",sep="",paste("_",NrList,sep="",paste(".csv",sep="")))
TogetherName <- rep(TogetherName,each=2)

Outcome <- "OutcomeMod"
OutcomeAllFin <- cbind(TogetherName,Outcome,OutcomeAll)

######################################################################################################

#Save together
save.xlsx(paste("Results/MoransITable_CoreModels",".xlsx",sep=""),PropensityAllFin,OutcomeAllFin)

#######################################################################################################
######################################################################################################

#2. Test for Endogeneity

##################################################################################################################
##################################################################################################################

#Endogeneity test for one model
test <- ResidCorreEndogeneityWElection("ProgrammeOutcome/ZH/Kcal/","1.ZHandKcal_core.R",
                               "#start2_","#fin2_","2004","2004to2013","ZH",
                               sqrtFuncEmpty,"ProgrammeOutcome/ZH/Kcal/",
                               "AllCoreMods_W_Election_Full.rda",NA,"Biomefor2004Analysis",
                               "AllCoreMods_W_Election_NoInfl.rda",3)


#Endogeneity test for all core models
sourcePartial("ZH_RobustnessTests.R",startTag="#starttest",endTag="#endtest")

CorrResEndog1 <- mapply(function(x,y,z,w,a,b,c,d,e,f,g,h) ResidCorreEndogeneityWElection(x,y,z,w,a,b,c,d,e,"AllCoreMods_W_Election_Full.rda",
                                                                                         f,g,"AllCoreMods_W_Election_NoInfl.rda",h),Folderlist,Scriptlist,
                        startList,endList,
                        YearSampList,TimeperiodList,ProgList,sqrtL,FolderSavelist,nrL,biomeL,NrList,SIMPLIFY = FALSE)

CorrResEndog1tab <- data.frame(do.call("rbind",CorrResEndog1))

##################################################################################################################

#Save
save.xlsx(paste("Result/Endogeneity_CoreModels",".xlsx",sep=""),CorrResEndog1tab)

##################################################################################################################
